全球变化情景下全球水稻高温致灾风险模拟栅格
数据集（2030s, 2050s）

苏  鹏1,2，王静爱1,3*，张安宇1，王  然1

1. 北京师范大学地理科学学部，北京 100875;
2. 青海师范大学地理科学学院，西宁 810016；
3. 青藏高原地表过程与生态保育教育部重点实验室，西宁 810008；

摘  要：本文利用MaxEnt模型、供需重分配方法和历史时期（1970–2000年）水稻种植区及产量数据，对未来近期2030s（2016–2035年）和未来中期2050s（2046–2065年）不同升温情景下（rcp2.6–ssp1，rcp4.5–ssp2和rcp8.5–ssp3）水稻的潜在适宜区进行模拟；然后，与计算得到的未来不同升温情景下的极端高温致灾叠加，得到水稻极端高温暴露区域数据。通过EPIC模型计算水稻面向高温胁迫的脆弱性曲线，并通过未来升温情景下的高温致灾计算，以此作为水稻高温产量损失数据，构成全球变化情景下全球水稻高温致灾风险模拟栅格数据集（2030s, 2050s）。该数据集内容包括：（1）历史时期的水稻种植区与3种升温情景下模拟得到的未来近期和未来中期的水稻潜在适宜种植区；（2）历史时期受极端高温影响的水稻分布区（即水稻高温暴露区）与3种升温情景下模拟得到的水稻高温暴露区；（3）历史时期极端高温导致的水稻产量损失（即水稻高温脆弱性）与3种升温情景下模拟得到的水稻产量损失。水稻种植区与高温暴露区数据空间分辨率为0.25°×0.25°，脆弱性数据空间分辨率为0.5°×0.5°。数据集存储为.tif格式，由21个数据文件组成，数据量为46.8 MB（压缩为1个文件，3.67 MB）。

关键词：全球变化；水稻；作物分布；暴露；脆弱性；产量损失

DOI: https://doi.org/10.3974/geodp.2022.04.05

CSTR: https://cstr.escience.org.cn/CSTR:20146.14.2022.04.05

数据可用性声明：

本文关联实体数据集已在《全球变化数据仓储电子杂志（中英文）》出版，可获取：

https://doi.org/10.3974/geodb.2022.06.04.V1或https://cstr.escience.org.cn/CSTR:20146.11.2022.06.04.V1.

1  前言

最新IPCC特别报告指出，按照目前全球的变暖速率，到2040年全球气温将上升1.5 ℃[1]，这会使得未来水稻作物的生境条件将发生变化。在不考虑作物品种改良的情况下，基于对自然环境的适应，水稻在未来的潜在分布区将会发生变化。同样的，未来全球面向水稻的高温致灾也将发生空间位置上的改变，这使得需要一套面向未来水稻高温致灾的承灾体以及暴露与产量损失数据。从而能够直观反映水稻在未来不同情景与时期的种植分布，灾害损失情况，从而有助于因地制宜地制定措施减少损失[2]，对世界粮食安全具有重要意义。

2  数据集元数据简介

《全球变化情景下全球水稻高温致灾风险数据集（2030s, 2050s）》[3]的名称、作者、地理区域、数据年代、时间分辨率、空间分辨率、数据集组成、数据出版与共享服务平台、数据共享政策等信息见表1。

 

表1  《全球变化情景下全球水稻高温致灾风险数据集（2030s, 2050s）》元数据简表

条  目

描  述

数据集名称

全球变化情景下全球水稻高温致灾风险模拟栅格数据集（2030s, 2050s）

数据集短名

GlobalRiceRisk

作者信息

苏鹏ABH-3399-2021, 北京师范大学地理科学学部，青海师范大学地理科学学院, 201947331031@stu.qhnu.edu.cn

王静爱AAA-6406-2022, 北京师范大学地理科学学部，青藏高原地表过程与生态保育教育部重点实验室，jwang@bnu.edu.cn

张安宇AAA-6787-2022, 北京师范大学地理科学学部, zay@mail.bnu.edu.cn

王然AAE-1120-2019, 北京师范大学地理科学学部, wangr0225@163.com

地理区域

全球

数据年代

2000s、2030s和2050s

时间分辨率

1970–2000（2000s）、2016–2035（2030s）、2046–2065（2050s）的时段内多年平均值

空间分辨率

0.25°Í0.25°，其中产量损失数据分辨率为0.5°Í0.5°

数据格式

.tif

数据量

46.8 MB

数据集组成

未来水稻种植面积、极端高温下的水稻暴露面积、极端高温下的水稻产量损失。

基金项目

中华人民共和国科学技术部（2016YFA0602402）

数据计算环境

ArcGIS、Matlab

出版与共享服务平台

全球变化科学研究数据出版系统 http://www.geodoi.ac.cn

地址

北京市朝阳区大屯路甲11号100101，中国科学院地理科学与资源研究所

数据共享政策

全球变化科学研究数据出版系统的“数据”包括元数据（中英文）、通过《全球变化数据仓储电子杂志（中英文）》发表的实体数据集和通过《全球变化数据学报（中英文）》发表的数据论文。其共享政策如下：（1）“数据”以最便利的方式通过互联网系统免费向全社会开放，用户免费浏览、免费下载；（2）最终用户使用“数据”需要按照引用格式在参考文献或适当的位置标注数据来源；（3）增值服务用户或以任何形式散发和传播（包括通过计算机服务器）“数据”的用户需要与《全球变化数据学报（中英文）》编辑部签署书面协议，获得许可；（4）摘取“数据”中的部分记录创作新数据的作者需要遵循10%引用原则，即从本数据集中摘取的数据记录少于新数据集总记录量的10%，同时需要对摘取的数据记录标注数据来源[4]

数据和论文检索系统

DOI，CSTR，Crossref，DCI，CSCD，CNKI，SciEngine，WDS/ISC，GEOSS

3  数据研发方法

3.1  基础数据

基础数据是本研究所用的原始数据，包括孕灾环境数据（全球海拔和土壤理化性质），承灾体数据（水稻种植范围、面积）和情景数据（气候数据）。本文选取USGS的全球数字高程模型以及GAEZ的全球地表坡度数据作为地形数据；全球土壤数据是来自ISRIC的全球土壤理化性质数据；水稻收获面积数据来自EARTHSTAT、SPAM 2005 v2.0、MIRCA2000；气候数据是来自NASA的全球地表交换降尺度逐日预测，包括降水、最高温与最低温，且使用多模式数据结果进行平均得到多模式综合结果。详细数据列表见表2。

 

表2  基础数据信息及来源统计表

类型

名称

年份

规格

来源

孕灾

环境

Global Multi-resolution Terrain Elevation Data (GMTED2010)

2010

1 km×1 km

United States Geological Survey (USGS)
https://topotools.cr.usgs.gov/gmted_viewer/

WISE derived soil properties on a 30 by 30 arc-seconds global grid

2012

30″×30″

International Soil Reference and Information Centre (ISRIC)
http://www.isric.org

承灾体

水稻种植范围、面积和产量数据：

①Harvested Area and Yield for 175 Crops year 2000

②SPAM 2005 v2.0

③MIRCA2000

2000或2005

5′×5′

①http://www.earthstat.org/harvested-area-yield-
175-crops/

②http://mapspam.info/maps/

③http://www.uni-frankfurt.de/45218031/data_
download?

FAO水稻产量、收获面积和单产统计数据

1960–2015

国家单元&亚国家单元

http://www.fao.org/faostat/en/#data/QC；http://kids.fao.org/agromaps/

FAO水稻需求量（来源于食物平衡表）

1960–2015

国家单元

http://www.fao.org/faostat/en/#data/FBS

FAO人口统计数据

1960–2015

国家单元

http://www.fao.org/faostat/en/#data/OA

FAO GDP统计数据

1960–2015

国家单元

http://www.fao.org/faostat/en/?#data/MK

水稻生育期

1961–1990

5’×5’

https://nelson.wisc.edu/sage/data-and-models/crop-calendar-dataset/index.php

情景

RCPs：NASA全球地表交换降尺度逐日预测（最高温、最低温、降水）

1960–2099

0.25°×0.25°

https://cds.nccs.nasa.gov/nex-gddp/

SSPs：全球格网化人口和GDP情景数据集

2005–2099

0.25°×0.25°

http://www.cger.nies.go.jp/gcp/population-and-gdp.html, Murakami, D. and Yamagata, Y. (2016)

Land-Use Harmonization2（LUH2）

1970–2100

0.25°×0.25°

http://luh.umd.edu/data.shtml

WorldClim V1.4历史风速

1970–2000

5′×5′

http://www.worldclim.com/version1

未来情景风速（4个模式）

1970–2099

0.5°×0.5°

https://www.isimip.org/gettingstarted/input-data-bias-correction/

 

3.2  算法原理

3.2.1  水稻种植分布预估

3.2.1.1  水稻潜在适宜区预估

本数据使用MaxEnt模型预估未来种植适宜性预估使用，该模型基于最大熵方法，对物种生态位和分布进行建模。模型的输入数据分为作物分布样本以及作物生境变量，参考本研究之前的工作[5]，本文的样本点数据及生境变量选择如下。

首先，本文将Harvested Area and Yield for 175 Crops year 2000数据与各栅格面积相除，得到作物种植面积占栅格面积的比率数据。其次，将作物种植面积占栅格面积的比率划分出四组，按适宜度最高组中的种植面积与抽样点比例在其他三组中确定点位选取的数量，并利用MATLAB中的randperm函数对各组样本点进行随机排序，按排序从前到后选择规定数量的样本点。选取样本点总数为2,228，四组样本点选择数量如下：559（面积<10%），768（10%≤面积<33%），541（33%≤面积<66%）以及360（66%≤面积）。

为了减少随机选点对结果的影响，本文随机选择了30套样本点。利用MaxEnt模型，采用相同的设置和生境变量，对这30套样本点分别计算适宜区结果后，得到30套全球水稻适宜区结果图，将30套结果数据取平均值，得到最终的适宜性结果。

本研究认为地形、气候[6–8]和土壤[9–11]是影响水稻自然适宜性最重要的三大因素。故选择高程、气候指标（22种）和土壤指标（20种）作为有可能对水稻自然适宜性产生影响的候选指标[12]，为了提取出对水稻的生长范围有明显影响的指标，本文对指标进行三次筛选。筛选后用于模型的指标见表4。

 

表3  MaxEnt模型输入指标统计表

数据类别

变量名称

变量含义

气候指标

BIO1

年均温

BIO2

平均日较差

BIO3

等温性

BIO5

暖月最高温

BIO8

湿季均温

BIO12

年降水量

BIO18

暖季降水量

Solar radiation

太阳辐射

Wind speed

风速

土壤指标

CECS

阳离子交换容量

CFRAG

土壤粗颗粒百分比

CNrt

碳氮比

GYPS

石膏含量

ORGC

有机碳含量

TAWC

体积含水量

TEB

盐基交换总量

地形指标

DEM

全球数字高程模型

 

3.2.1.2  水稻产量再分配预估

（1）水稻产量预估

水稻的需求量通过建立历史时期人均食用水稻供给。因其实际反映的是当年该国的水稻食物需求，故后称人均食用水稻需求与人均GDP之间的关系，来推算未来各SSP情景下水稻的消耗量。拟合函数参考前人的研究结果主要有幂函数线性模型[13]和对数线性模型[14]（拟合样例见图1）。得到未来人均食用水稻需求后，由于各国将水稻作为食物的比例不同，需要除以人均食用水稻需求占总水稻的比例以获得总水稻的需求量。这里首先定义水稻的食物需求量与国家的水稻总消耗量的比例不随时间发生变化，故该比例可由历史时期多年平均比例计算得到。最后将人均水稻需求量与人口相乘，即可得到未来各国的水稻需求量。

 

图1  人均GDP–人均水稻年需求量曲线拟合示例图

 

（2）水稻产量分配

水稻产量再分配的总体思路是根据各国水稻供给和需求量，使供需达到平衡。供给方面，使用适宜性变化预测出初始水稻收获面积，作为迭代的起始值，通过比较国家水稻的供给和需求，对水稻收获面积进行调整，如果供大于求，则减少水稻收获面积，反之则增加收获面积，直到供需平衡或所有耕地种植水稻为止。该模型可以表达为公式1–4同时成立的多目标优化模型。

                                                                                                                                (1)

式中，D，S，I分别代表国家水稻需求，供给和净进口量；

                                                                                                                           (2)

式中，，分别代表第i个栅格的水稻收获面积和C3作物收获面积。

                                                                         (3)

式中，n为迭代次数，n=1为根据适宜区变化推断的初始情况，Su为适宜性。

                                                                                                                    (4)

式中，为FAO预测的国家水稻灌溉率，，分别为栅格i中水稻灌溉面积和水稻收获面积。

图2详细展示了供给测的流程，该流程中考虑了灌溉因素。在水稻收获面积调整时，遵循以下三个原则：①使国家内的水稻供需达到平衡；②使国家内水稻灌溉比例尽可能靠近FAO预测值；③在增加/减少收获面积时，适宜区被等间距分为100级，从高适宜区开始增加/低适宜区开始减少，直到水稻种满该适宜区内所有耕地/该适宜区范围内无水稻种植。在实际操作时，对于增加水稻收获面积的国家，由于无法确定特定栅格增加范围内的灌溉比例，假设所有水稻增加面积的灌溉比例相同。

图2  未来水稻收获面积、产量预估框架流程图

 

3.2.2  暴露及产量损失预估方法

3.2.2.1  致灾因子计算

水稻的极端高温暴露用高温发生日数和超过阈值积温衡量，具体为生育期内温度高于水稻适宜生长的最高温度[15]（水稻的极端高温阈值为38 ℃）的日数和积温。本研究将单日平均气温超过38 ℃（≥38 ℃）定义为一次水稻高温暴露事件，并将生育期内高温暴露事件的累积胁迫值（GHTS）作为高温致灾强度。

3.2.2.2  暴露计算

水稻的暴露范围为水稻分布范围和致灾范围在空间上的叠加，而范围中填充的数值则是暴露的量，本数据所使用的暴露量指标为年水稻收获面积和水稻年产量。对于极端高温下的暴露，暴露范围为极端高温日数>0的范围与水稻种植范围的交集构成，暴露范围内显示水稻的收获面积。

3.2.2.3  脆弱性计算

脆弱性曲线是通过建立致灾因子强度与承灾体损失的二元关系表达脆弱性概念。本研究采用EPIC模型分别计算高温致灾指数与产量损失率，构建水稻高温脆弱性曲线。

EPIC模型是综合气候、土壤、水分、田间管理等因素的动力学模型，常用于作物产量模拟。本文运用EPIC模型计算得到高温致灾强度指数与作物损失率间的脆弱性曲线。计算公式如下。

本研究中的高温致灾强度指数（HSI）定义为：

                                                                                                                     (5)

式中，TSi为第i日温度胁迫值，n为生育期日数，TS为某情景下的生育期温度累积胁迫，TSmax为生育期累积胁迫的潜在最大值。HSI数值越大代表高温致灾强度越高。

本研究中的产量损失率（YL）定义为

                                                                                                             (6)

式中，Y为某情景下的模拟产量，Ymax为最优情景（TS=0）下的模拟产量。

根据上述计算得到HSI与相对应的YL通过logistics回归方程拟合脆弱性曲线，本研究中所采用的拟合方程为

                                                                     (7)

式中，a、b、c、d为曲线函数参数。

3.2.2.4  产量损失计算

在暴露区域内采用信息扩散理论计算得到高温致灾强度指数的概率密度曲线，将其与对应的脆弱性曲线相结合，得到水稻产量损失率的概率密度曲线。计算其损失期望，反映某时期内高温致灾影响下水稻产量损失的平均状态。

通过EPIC模拟某时期内各评价单元内的温度累积胁迫值，采用式5计算得到逐年高温致灾强度指数（HSI）。以此数据作为信息扩散样本，采用正态扩散方式实现概率密度分布估计。具体计算过程如下：

令U={u1,u2,…,un}是包含HSI可能值的离散论域，HSI的取值范围为0-1，本研究中论域的分辨率为0.000,1，因此U={0, 0.0001, 0.0002,…,1}。通过信息扩散函数（式8）将每个网格内的HSI所携带的信息扩散到每个ui中。

                                                                                  (8)

式中，k是每个网格的编码，h为正态扩散系数，可以通过式9计算得到。

                                                                                                 (9)

式中，a，b分别为HSI的最小值和最大值，m为样本数量（即年份）。则样本的信息累积和正态信息分布可以分别用式10与式11计算。

                                                                                                                     (10)

                                                                                                           (11)

式中，Ck是第k各样本的信息积累，F(HSIk,uj)为HSI的标准化信息分布，对于每个点uj，将所有归一化信息相加，得到来自给定样本HSI的uj处的信息增益。信息增益显示在式12中。

                                                                                                        (12)

通过对q(ui)的求和得到样本的扩散信息量：

                                                                                                                         (13)

进而计算HSI的概率密度分布：

                                                                                                                       (14)

定义为某段时期内水稻高温产量损失率的期望值（式15）。

                                                                                                               (15)

式中，p(uj)为高温致灾强度指数为uj时的概率估计值。

3.3  技术路线

本文数据分为三个部分：未来水稻种植面积数据，极端高温下水稻面积暴露数据，极端高温下水稻产量损失数据。其中未来水稻种植面积数据由自然与社会经济双情景计算得到，极端高温下水稻面积暴露数据由致灾数据与水稻种植分布数据共同计算得到，极端高温下水稻产量损失数据由EPIC模型模拟得到的水稻高温脆弱性曲线与致灾数据计算得到。详细流程见图3。

 

图3  数据计算技术路线框架图

4  数据结果与验证

4.1  数据集组成

数据在时间上分为1970–2000（历史时期），2016–2035（未来近期），2046–2065（未来中期）。模式分为rcp2.6–ssp1，rcp4.5–ssp2和rcp8.5–ssp3三种组合。

数据分辨率：0.25°×0.25°，其中极端高温下水稻产量损失的数据空间分辨率为0.5°×0.5°。

数据格式：.tif。

数据命名方式：类型_作物_时期_情景.tif。

4.2  数据结果

4.2.1  未来水稻种植分布数据

本文通过MaxEnt模型模拟水稻的自然适宜区，并在此基础上进行国家单元的供需分配，最终得到水稻在不同情景下的种植分布，并进行制图（图4）。通过两个时期的对比，相较于历史时期，在RCP8.5–SSP3情景下未来中期水稻的分布面积总体上呈现减少趋势，南美洲与印度半岛种植面积减少最为明显。

4.2.2  未来水稻极端高温暴露数据

将水稻种植分布范围与极端高温致灾范围进行叠加，得到水稻在各情景下极端高温暴露情况，并进行制图（图5）。通过两个时期的对比，相较于历史时期，在RCP8.5–SSP3情景下未来中期水稻的极端高温暴露面积总体上呈现增加趋势，非洲与我国南方地区的水稻极端高温暴露面增加最为明显。

4.2.3  未来水稻面向极端高温产量损失数据

将水稻面向高温致灾的脆弱性曲线与未来情景下的高温致灾结合，得到水稻在各情景下受极端高温致灾造成的产量损失比率，并进行制图（图6）。通过两个时期的对比，相较于历史时期，在RCP8.5–SSP3情景下未来中期极端高温致灾对水稻造成损失在分布以及在程度上均有增大。

4.3  数据结果验证

本数据集基于致灾（水稻高温事件）以及承灾体（水稻种植分布）计算得到，其中致灾数据由CMIP5的各气候模式计算得到，这些气候模式数据的可用性已得到广泛地证明[16]。水稻种植分布计算分为自然适宜区预估以及水稻种植再分配，其中种植再分配后的水稻种植区是基于历史时期水稻种植数据模拟得到，因此难以使用历史时期水稻收获数据进行验证。因此本数据集精度验证主要集中于水稻自然适宜区的验证。在之前的研究中，我们通过遥感影像分类等方法对适宜区进行了验证[6]，证明了数据的合理性。在此我们将对样本点的随机选择进行更深入的讨论。

因为当随机样本组数量大于30时，波动较大区域的不确定性降低幅度明显减弱。因此综合考虑降低不确定性和效率因素，我们选择将随机样本组的数量定为30组，用以降低结果的不确定性。本文逐网格计算30组随机样本组的标准差，并将其进行空间制图（图7）。

图4  历史时期与综合气候模式下RCP8.5–SSP3情景下未来中期水稻种植分布图

 

图5  综合气候模式下的历史时期与RCP8.5–SSP3情景下未来中期水稻极端高温暴露分布图

 

图6  综合气候模式下的历史时期与RCP8.5–SSP3情景下未来中期水稻极端高温产量损失分布图

 

将标准差进行空间化制图，发现标准差较大区域主要集中在帕米尔高原以及小亚细亚半岛，表明这些区域水稻分布模拟结果在不同样本组间差异较大，而这两个区域内的实际水稻种植面积均较小，因此对整体模拟结果影响不大。

5  讨论和总结

本数据集在利用MaxEnt模型得到不同升温情景下的水稻适宜区，并在此基础上进行社会经济情景的水稻种植再分配，得到未来不同升温情景下的水稻种植分布，并以此作为

图7  样本组（共计30组）标准差分布图

 

研究的承灾体。通过遥感影像分类、模型ROC曲线、统计验证、模型模拟预实验等方法，本文验证了承灾体预估结果的有效性。通过计算未来极端高温致灾，并将其与承载体叠加，得到未来极端高温致灾下的水稻暴露面积数据。且利用EPIC模型进行模拟得到不同强度高温胁迫下的水稻产量损失率，以此得到水稻面向高温致灾的脆弱性曲线。在此基础上，得到未来不同时期相较于最优产量的水稻产量损失率。

本数据集通过综合考虑自然环境与社会经济双重影响，预估得到未来水稻种植分布。并在此基础上进行承灾体与致灾因子双动态变化综合考虑的暴露与产量损失数据计算。本数据集一方面对全球气候变化背景下未来水稻种植、粮食安全研究有显著意义，另一方面，对气候变化影响评估也有重要价值，对于农业领域的气候变化影响评估所需的水稻空间分布资料，本数据集也可提供重要参考。

 

作者分工：王静爱和苏鹏对数据集的开发做了总体设计；张安宇和王然采集和处理了原始数据；苏鹏、张安宇和王然设计了模型和算法；苏鹏做了数据验证。

 

利益冲突声明：本研究不存在研究者以及与公开研究成果有关的利益冲突。

 

参考文献

[1]      Intergovernmental Panel on Climate Change. Global warming of 1.5 ℃ [R]. Geneva, Switzerland: World Meteorological Organization, 2018.

[2]      Fraser, S., Wood N., Johnston D., et al. Variable population exposure and distributed travel speeds in least-cost tsunami evacuation modelling [J]. Natural Hazards and Earth System Sciences, 2014, 2(6): 4163–4200.

[3]      苏鹏, 王静爱, 张安宇等. 全球变化情景下全球水稻高温致灾风险模拟栅格数据集(2030s, 2050s)[J/DB/OL]. 全球变化数据仓储电子杂志, 2022. https://doi.org/10.3974/geodb.2022.06.04.V1. https://cstr.escience.org.cn/CSTR:20146.11.2022.06.04.V1.

[4]      全球变化科学研究数据出版系统. 全球变化科学研究数据共享政策[OL]. https://doi.org/ 10.3974/dp.policy.2014.05 (2017年更新).

[5]      Su, P., Zhang, A., Wang, R., et al. Prediction of future natural suitable areas for rice under representative concentration pathways (RCPs) [J]. Sustainability, 2021, 13(3): 1580.

[6]      Yoshida, S., Parao, F. T., Climatic influence on yield and yield components of lowland rice in the tropics [J]. Climate and Rice, 1976, 20: 471–494.

[7]      Lobell, D. B., Schlenker, W., Costa-Roberts, J. Climate trends and global crop production since 1980 [J]. Science, 2011, 333(6042): 616–620.

[8]      Lobell, D. B., Gourdji, S. M. The influence of climate change on global crop productivity [J]. Plant physiology, 2012, 160(4): 1686–1697.

[9]      Azam, F., Comparative effects of organic and inorganic nitrogen sources applied to a flooded soil on rice yield and availability of N [J]. Plant & Soil, 1990, 125(2): 255–262.

[10]   Kravchenko, A. N., Bullock, D. G. Correlation of corn and soybean grain yield with topography and soil properties [J]. Agronomy Journal, 2000, 92(1): 75–83.

[11]   Jiang, X., Xie, D. Combining ridge with no-tillage in lowland rice-based cropping system: long-term effect on soil and rice yield [J]. Pedosphere, 2009, 19(4): 515–522.

[12]   Wang, R., Jiang, Y., Su, P., et al. Global spatial distributions of and trends in rice exposure to high temperature [J]. Sustainability, 2019, 11: 6271.

[13]   Tilman, D., Balzer, C., Hill, J., et al. Global food demand and the sustainable intensification of agriculture [J]. Proceedings of the National Academy of Sciences of the United States of America, 2011, 108(50): 20260–20264.

[14]   Bijl, D. L., Bogaart, P. W., Dekker, S. C., et al. A physically-based model of long-term food demand [J]. Global Environmental Change, 2017, 45: 47–62.

[15]   Watson, R. T., Zinyowera, M. C., Moss, R. H., Impacts, adaptations and mitigation of climate change: scientific-technical analyses [R]. Contribution of Working Group II to the Second Assessment Report of the Intergovernmental Panel on Climate Change, University Press, 1996.

[16]   Taylor, K. E., Stouffer, R. J., Meehl, G. A. an overview of CMIP5 and the experiment design [J]. Bulletin of the American Meteorological Society, 2012, 93(4): 485–498.